Taxonomic filtering

11 January 2022

Read the taxonomy files back into Rstudio after performing the blastn search.

library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)

pool_taxonomy <-read.delim("blast_outputs/12S_poolseqs_taxlineage.txt", header = FALSE)

head(pool_taxonomy)

clean up the header a bit

# use the full taxonomy rather than the seq id to collapse identical entries
pool_tax_df <- pool_taxonomy %>%
  filter(V4 > 100) %>% # make sure all retained matches are >100 bp
  select(-V2, -V5, -V6, -V7, -V14, -V8, -V9, -V10, -V11, -V12, -V15, -V16) %>%  #remove unnecessary columns
  group_by(V1, V17) %>% # group by the sequence key and the full taxonomy to reduce duplicate entries
  unique() %>% # doing that reduced the number of entries from 146k to 17k
  rename(qseqid=V1, perc_id=V3, length=V4, sscinames=V13, taxonomy=V17) %>% #rename headers
  filter(!str_detect(sscinames, "environmental")) %>% # filter out any environmental samples
  filter(!str_detect(sscinames, "synthetic")) %>% # filter out any synthetic "samples"
  filter(perc_id >= 98) # seems like some of the matches below 98% are dubious (jellyfish and herring <1% different??)
 

Maybe break up that taxonomic info so that we can focus on the things that we can identify to some relevant level (at least order!?)

# first figure out where there are issues with that taxonomy format:
genus.hybrids <- pool_tax_df %>%
  filter(!str_detect(taxonomy, ";")) %>% 
  # looking at these, they're probably hybrids or not good species, which is why the taxonomy broke down. 
  separate(sscinames, into = c("genus"), extra = "drop", remove = FALSE)

# and where there aren't issues
taxon_df <- pool_tax_df %>%
  filter(str_detect(taxonomy, ";")) %>%
  separate(taxonomy, into=c("kingdom", "phylum", "class", "order", "family", "genus", "species"), sep = ";")

# create a separate taxonomic lineage df with unique entries
lineage_df <- taxon_df %>%
  ungroup() %>%
  select(kingdom, phylum, class, order, family, genus) %>%
  unique()


# there are only 18 unique genera in this small messed up part of the dataset 
genus_fixed <- genus.hybrids %>%
  ungroup() %>%
  select(genus) %>%
  unique() %>% # can I fix the lineage by adding existing data?
  left_join(lineage_df, by = "genus") %>% # now the Tachyeres duck is the only loner
# fix that manually
  mutate(family = ifelse(genus == "Tachyeres", 
                         "Anatidae", family)) %>%
    mutate(order = ifelse(genus == "Tachyeres",
                        "Anseriformes", order)) %>%
      mutate(class = ifelse(genus == "Tachyeres",
                        "Aves", class)) %>%
        mutate(phylum = ifelse(genus == "Tachyeres", 
                               "Chordata", phylum)) %>%
          mutate(kingdom = ifelse(genus == "Tachyeres",
                                    "Eukaryota", kingdom))

# add the lineage info back to the rest of the data
genus_fixed_df <- genus.hybrids %>%
  ungroup() %>%
  left_join(genus_fixed, by = "genus") %>%
  select(-taxonomy)

# and merge them together
fixed_tax_df <- taxon_df %>%
  bind_rows(genus_fixed_df) %>%
  # and fix the species syntax with underscores
   mutate_if(is.character, str_replace_all, ' ', '_') %>%
   filter(!str_detect(sscinames, "_x_"))

`mutate_if()` ignored the following grouping variables:
Column `qseqid`

  # remove hybrids because the mtDNA sequence actually is only one or the other of the parent species

Sort and clean-up the results based on % identity

# how many total sequences?
fixed_tax_df %>%
  select(qseqid) %>%
  unique()

NA

There are four categories: 1. sequences that match a single species unambiguously (the minority)

Sequences that match multiple species are divided in three categories: 2. top matches > 2% identity than second-ranked matches 3. top matches < 2% identity than second-ranked matches 4. Multiple top matches with the same % identity

# 1. sequences that are unambiguously a single species
single_spp_seqs <- fixed_tax_df %>% 
  group_by(qseqid) %>%
  add_tally(name = "n_taxa") %>%
  filter(n_taxa == 1)
  

140 sequences are single-species

# remove the single-species seqs from the dataframe and then rank the hits by % identity for the remaining seqs
seq_id_diff <- fixed_tax_df %>%
  anti_join(., single_spp_seqs) %>%
  select(-length) %>%
  group_by(qseqid, species, genus, family, order, class, phylum, kingdom) %>%
    mutate(seq_percID = max(perc_id)) %>%
    group_by(qseqid, species, genus, family, order, class, phylum, kingdom, seq_percID) %>%
  summarise(max(seq_percID)) %>% # take just the top hit for each taxon (for each sequence)
  select(-`max(seq_percID)`) %>%
  ungroup() %>%
  group_by(qseqid) %>%
      mutate(id_rank = rank(desc(seq_percID), ties.method = "min")) %>% # rank the taxonomic hits per sequence by % id
       mutate(top_perc = max(seq_percID)) %>% # designate the highest % id for the best taxonomic hit in each sequence (in some, but not all cases, this is 100%)   
      mutate(diff = top_perc - seq_percID) %>% # calculate the difference between the % identity of the top hit and each subsequent taxonomic hit
        arrange(diff)

Joining, by = c("qseqid", "perc_id", "length", "sscinames", "kingdom", "phylum", "class", "order", "family", "genus", "species")
`summarise()` has grouped output by 'qseqid', 'species', 'genus', 'family', 'order', 'class', 'phylum', 'kingdom'. You can override using the `.groups` argument.

seq_id_diff %>%
  filter(diff > 0)

NA

Now I have the single best entry for each species for each sequence ranked and with the difference between the first and second ranked entries calculated.

For sequences with multiple top hits, where the difference between ranked taxa = 0, I will end up defaulting to genus- or family-level ID (or carrying the individual species info around in some capacity). I will do the same for any sequences where the difference betweeen the first and second ranked taxa is < 2%.

Figure out which differences are > 2% and eliminate those first?

# filter out any taxa that are >2% less matching identity than the top taxonomic hit for a given sequence
to_remove_low_perc_hits <- seq_id_diff %>%
  ungroup() %>%
  group_by(qseqid) %>%
  filter(diff > 2)

keepers <- seq_id_diff %>%
  anti_join(to_remove_low_perc_hits)

Joining, by = c("qseqid", "species", "genus", "family", "order", "class", "phylum", "kingdom", "seq_percID", "id_rank", "top_perc", "diff")

# this data frame includes only those taxonomic hits that should be considered.
# so now I need to determine whether they should be assigned to genus, family, order, etc. 
singletons <- keepers %>%
  select(qseqid) %>%
  tally() %>%
  filter(n == 1)

# these are the seqs that now have only a single match
singleton_df <- singletons %>%
  left_join(keepers) %>%
  select(-n) %>%
  bind_rows(single_spp_seqs) # combine the single spp data

## Genus-level matches
# remove the singletons from the bigger df 
single_genus <- keepers %>%
  anti_join(singleton_df)%>% # at best, these should be genus-level matches
  group_by(qseqid, genus) %>%
  tally() %>%
  ungroup() %>%
  group_by(qseqid) %>%
  tally() %>%
  filter(n == 1) %>% # seqs that match a single genus
  select(-n) %>%
  left_join(., keepers) %>%
  mutate(taxonomic_level = "genus") %>%
  mutate(taxon = genus)
  

## Family-level matches
single_family <- keepers %>%
  anti_join(singleton_df)%>%
  anti_join(single_genus) %>%
  group_by(qseqid, family) %>%
  tally() %>%
  ungroup() %>%
  group_by(qseqid) %>%
  tally() %>%
  filter(n == 1) %>% # seqs that match a single family
  select(-n) %>%
  left_join(., keepers) %>%
  mutate(taxonomic_level = "family") %>%
  mutate(taxon = family)


## Order-level matches
single_order <- keepers %>%
  anti_join(singleton_df)%>%
  anti_join(single_genus) %>%
  anti_join(single_family) %>%
  group_by(qseqid, order) %>%
  tally() %>%
  ungroup() %>%
  group_by(qseqid) %>%
  tally() %>%
  filter(n == 1) %>% # seqs that match a single order
  select(-n) %>%
  left_join(., keepers) %>%
  mutate(taxonomic_level = "order") %>%
  mutate(taxon = order)


## Class-level matches
single_class <- keepers %>%
  anti_join(singleton_df)%>%
  anti_join(single_genus) %>%
  anti_join(single_family) %>%
  anti_join(single_order) %>%
  group_by(qseqid, class) %>%
  tally() %>%
  ungroup() %>%
  group_by(qseqid) %>%
  tally() %>% 
  filter(n == 1) %>% # seqs that match a single class
  select(-n) %>%
  left_join(., keepers) %>%
  mutate(taxonomic_level = "class") %>%
  mutate(taxon = class)


## Phylum-level matches
single_phylum <- keepers %>%
  anti_join(singleton_df)%>%
  anti_join(single_genus) %>%
  anti_join(single_family) %>%
  anti_join(single_order) %>%
  anti_join(single_class) %>%
  group_by(qseqid, phylum) %>%
  tally() %>%
  ungroup() %>%
  group_by(qseqid) %>%
  tally() %>% 
  filter(n == 1) %>% # seqs that match a single phylum
  select(-n) %>%
  left_join(., keepers) %>%
  mutate(taxonomic_level = "phylum") %>%
  mutate(taxon = phylum)

# Kingdom level matches
single_kingdom <- keepers %>%
  anti_join(singleton_df)%>%
  anti_join(single_genus) %>%
  anti_join(single_family) %>%
  anti_join(single_order) %>%
  anti_join(single_class) %>%
  group_by(qseqid, kingdom) %>%
  tally() %>%
  ungroup() %>%
  group_by(qseqid) %>%
  tally() %>% 
  filter(n == 1) %>% # seqs that match a single kingdom
  select(-n) %>%
  left_join(., keepers) %>%
  mutate(taxonomic_level = "kingdom") %>%
  mutate(taxon = kingdom)

Seems… unlikely that a jellyfish would be 1% dissimilar from a herring. BUT, it’s below 98%, so maybe that’s at least one reason to be more stringent on the upstream side of this?

Adding that 98% filter in back up top.

Modify the singleton_df to include the right variable headers

single_spp <- singleton_df %>%
  select(-perc_id, -length, -sscinames, -n_taxa) %>%
  mutate(taxonomic_level = "species") %>%
  mutate(taxon = species)

Okay, so that df should contain the basic taxa in the dataset assigned to the taxonomic level for which we can be relatively confident (>2% difference between the top-and-next highest hit).

Filter the eDNA taxonomy data frame by Fish Atlas

I fundamentally have an issue with this step: we want to use biogeographic information to remove illogical matches, but we know that Fish Atlas is not comprehensive and so this “filter” is a very manual and tedious check. For example, sablefish, which is very common in nearshore SE AK, is not included in the Fish Atlas data.

Given that, here is how I will proceed:

Using family as the level of matching: First, I want to fill out the order for the fish atlas data using my data frame, then I’ll filter the data by family, but retain the species and genus designations as appropriate base on the prior taxonomic filtering step.

# read in fish atlas data from Pat:
fishatlas <- read.csv("../data/FAdata.csv") %>%
  # remove unnecessary columns
  select(Sp_CommonName, Sp_ScientificName, Fam_ScientificName) %>%
  arrange(Fam_ScientificName)

# clean-up the format (remove spaces)
fishb <- fishatlas %>% 
  mutate_all(~(str_replace(., " ", "_"))) %>%
  rename(common = Sp_CommonName, family = Fam_ScientificName)

# fix rockfishes so that it matches our eDNA data frame
fish.atlas.fixed <- fishb %>%
  mutate(family = ifelse(str_detect(common, "rockfish"), "Sebastidae", family))

quick check: which taxa are filtered by the fish atlas?

# which taxa are filtered out when using the fish atlas?
# it turns out some of these are definitely real detections.
just.fish %>%
  anti_join(., fish.atlas.fixed) %>%
  unique() %>%
  arrange(family) %>%
  filter(!str_detect(species, "Carassius")) %>%
  filter(!str_detect(species, "Cyprinus")) %>%
  #select(family) %>%
  unique() %>%
  arrange(qseqid) %>%
  filter(taxonomic_level %in% c("species", "genus", "family", "order")) %>%
  select(qseqid, species, family, taxon, taxonomic_level)

Joining, by = "family"

**that is where sablefish and rockfish are filtered out! Why?? Manual checks on the taxa filtered out by Fish Atlas.

Focus on the “taxon” field here, since those are the data that will be carried forward The Merluccius should be filtered out because of geography. Lepomis_macrochirus = bluegill Labidesthes_sicculus = brook silverside Leuciscidae = freshwater fish (Cyprinidae) Ambloplites = freshwater fish Merluccius = Atlantic/Europe Mallotus_villosus = capelin (possible in AK??) Micropterus_dolomieu = smallmouth bass Etheostoma = freshwater Engraulidae = most likely from the fish food Lipolagus_ochotensis = black-eared smelt (deep water… but in the N Pacific) Catostomidae = freshwater suckers Percidae = fresh/brackish waters of N America Anoplopomatidae = sablefish!!! Scomber_scombrus = Atlantic mackerel (fish food, maybe) Echiichthys_vipera = Atlantic Micropterus = black basses Macrouridae = deep water (Grenadiers or rattails) Pimephales = freshwater Percina_peltata = freshwater/Eastern US Trachinus_draco = Atlantic Bagre_marinus = catfish Anarrhichthys_ocellatus = wolf eel

Here’s where the painful part comes in - we want to keep the following species that are otherwise filtered out by the Fish Atlas join:

“Anoplopomatidae” “Anarrhichthys_ocellatus” “Macrouridae” “Lipolagus_ochotensis” “Mallotus_villosus”

Probably the simplest way to do this is to keep the inner join with Fish Atlas and then add these back on.

fishatlas.edna.tax <- fish.atlas.fixed %>%
  inner_join(., just.fish, by = "family") %>%
  mutate(to.filter = ifelse(taxonomic_level == "species" & species != Sp_ScientificName, "yes", "no"))  # this provides us a way to filter mis-identified taxa to the species level based on fish atlas data - but these taxa are fine to retain at the genus/family level

Quick question: how many taxa get removed by because the family matches but the species does not?

# remove duplicates so that I can take a closer look at which species-level IDs would be removed.
fishatlas.edna.tax %>%
  filter(to.filter == "yes") %>%
  select(-qseqid) %>%
  select(taxon) %>%
  unique()

NA
NA

In total, 8 species would be removed from the dataset due to matching family but mismatching species with Fish Atlas. Let’s verify each one.

Oncorhynchus_tshawytscha - Chinook salmon… definitely present in the marine nearshore outside of Juneau. Probably shouldn’t remove. Leptocottus_armatus - Pacific staghorn sculpin - also quite likely to be present Artedius_fenestralis - Padded sculpin - distribution extends to AK.. unknownn. Atheresthes_evermanni - Kamchatka flounder - unlikely Sardina_pilchardus - European pilchard - unlikely (fish food) Opisthonema_oglinum - Atlantic thread herring - unlikely (fish food) Hexagrammos_decagrammus - kelp greenling (definitely possible) Gymnammodytes_semisquamatus - Smooth sandeel - distribution = Europe (maybe fish food)

So the conclusion is that 50% of those are possible or likely in our study system and we should keep the following: Oncorhynchus_tshawytscha Leptocottus_armatus Artedius_fenestralis Hexagrammos_decagrammus Atheresthes_evermanni -should be switched to genus-level because arrowtooth flounder is a sister-species of Kamchatka flounder

I hate doing this, but I’m going to have to manually filter the four to preserve the others:

# filter based on fish atlas, but retain species that were double-checked
fishatlas.edna.filtered <- fishatlas.edna.tax %>%
  filter(!taxon %in% c("Sardina_pilchardus", "Opisthonema_oglinum", "Gymnammodytes_semisquamatus")) %>%
  # and then add back on the taxa from extra.keepers
  bind_rows(extra.keepers) %>%
  # switch Atheresthes_evermanni to a genus-level match
  mutate(taxon = ifelse(species == "Atheresthes_evermanni", 
                        "Atheresthes", taxon)) %>%
  mutate(taxonomic_level = ifelse(species == "Atheresthes_evermanni",
                        "genus", taxonomic_level)) %>%
  select(-to.filter)

And just to confirm that that took care of the outstanding issues, look for sebastes, sablefish, etc.

# export that fish-only taxonomy df
fishatlas.edna.filtered %>%
  write.csv("csv_outputs/12Sfishatlas.edna.taxonomy.csv")

Now, looking at the sample/sequence dataframe

There are only 748 unique sequences that made it through the blast/taxonomy gauntlet, so why are there an order of magnitude more sequences in the seqtab.nochim data tables. Maybe because they didn’t make it through the BLAST filtering? i.e., no match at >98% identity?

1. I need to make the ASV-sample table usable without the actual sequences as the headers.
2. read in the new, clean metadata (that I pushed to the VM) and connect it with the sample info
3. bung in all together into a df and send to Wes for scrutiny.

STOP HERE

go to 06-merge-samples-w-taxonomy.Rmd instead for the latest version of this analysis.

Read in ASV-sample table for Poolseq data

asv_tbl <- read.csv("csv_outputs/12S_poolseq_ASVtable.csv") %>%
  rename(sample.PCR = X) %>%
  separate(sample.PCR, into = c("sample", "rep"), remove = F)

Combine taxonomy and sample ASV table

# pivot longer for ease of joining
asv_long_tbl <- asv_tbl %>%
  pivot_longer(cols = 4:7705, names_to = "ASV", values_to = "count")

# join data frames
uncollapsed_poolseq_tax_df <- asv_long_tbl %>%
  left_join(., sorted_tax_df, by = c("ASV" = "qseqid"))

# collapse these based on the appropriate taxonomic level
collapsed_poolseq_tax_df <- uncollapsed_poolseq_tax_df %>%
  select(sample.PCR, sample, rep, ASV, count, taxon, taxonomic_level) %>% ## I need to go back and fix the missing % for the species-level hits
  unique()

That data frame still contains a bunch of NAs that didn’t have blast hits. Let’s go ahead and remove them.

poolseq_tax_tidydf <- collapsed_poolseq_tax_df %>%
  filter(!is.na(taxon)) %>% # remove NAs
  filter(count > 0) # remove entries with zero reads for that sequence

Add metadata to that tidy frame

# metadata tidied and ready to go
meta <- read.csv("../metadata/amalga_clean_metadata.csv")

poolseq_tax_tidymeta <- poolseq_tax_tidydf %>%
  left_join(., meta, by = "sample") #%>%
  #filter(!is.na(label)) # remove aquarium samples?

poolseq_tax_filtered_w_meta <- poolseq_tax_tidymeta %>%
  arrange(distance) %>%
  filter(taxonomic_level != "kingdom") %>% # remove these totally uninformative matches
  # collapse multiple entries for the same taxon in the same sample replicate
  group_by(sample.PCR, taxon) %>%
  mutate(count = sum(count)) %>%
  select(-ASV) %>% 
  unique() # that reduces the number of rows from 7k to 4k

That’s pretty useful.

Quick check on negative controls

Let’s break it down in a few different ways:

# negatives?
poolseq_tax_filtered_w_meta %>%
  filter(label == "blank", taxon != "Homo_sapiens") %>%
  ggplot(aes(x=sample.PCR, y=count, fill=taxon)) +
    geom_bar(stat="identity") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90),
    axis.title.x = element_text(margin = margin(t=10))
  ) +
  labs(
    y = "reads",
    x = "extraction blanks"
  )

ggsave("pdf_outputs/12S_extraction_blanks.pdf")

Saving 7.26 x 4.49 in image

# pcr negatives
poolseq_tax_filtered_w_meta %>%
  filter(str_detect(sample.PCR, "NEGATIVE"), taxon != "Homo_sapiens") %>%
  ggplot(aes(x=sample.PCR, y=count, fill=taxon)) +
    geom_bar(stat="identity") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90),
    axis.title.x = element_text(margin = margin(t=10))
  ) +
  labs(
    y = "reads",
    x = "PCR negative controls"
  )

ggsave("pdf_outputs/pcr_negatives.pdf")

Saving 7.33 x 4.52 in image

# these PCR negatives look awesome.

# positives?
poolseq_tax_filtered_w_meta %>%
  filter(str_detect(sample.PCR, "POSITIVE"), taxon != "Homo_sapiens") %>%
  ggplot(aes(x=sample.PCR, y=count, fill=taxon)) +
    geom_bar(stat="identity") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90),
    axis.title.x = element_text(margin = margin(t=10))
  ) +
  labs(
    y = "reads",
    x = "PCR positive controls"
  )

I guess the goldfish doesn’t show up at any taxonomic level higher than class?

# species/genus/family-level matches?
poolseq_tax_filtered_w_meta %>%
  filter(!is.na(label), 
         taxon != "Homo_sapiens",
         taxon != "Canis_lupus",
         taxonomic_level %in% c("species")) %>%
  ggplot(aes(x=sample, y=count, fill=taxon)) +
    geom_bar(stat="identity") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90),
    axis.title.x = element_text(margin = margin(t=10))
  ) +
  labs(
    y = "reads",
    x = "sample replicates"
  )

NA
NA

Next up: organize by distance and look at composition. Also - could summarize/further collapse taxa by sample when sequences fall in multiple taxonomic levels.

stop here.

Run an NMDS?

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.5-7

To run an ordination you will need a data-frame consisting of plot by species (or trait) matrix AND a “groups” data-frame which should consist of plots with a coding variable for what group each plot belongs to, this will be used for plotting the ordination.

This is not really the format that we currently have, so let’s move things around a bit.

# we don't need the metadata for this - simplify
pool.vegan.df <- poolseq_tax_filtered_w_meta %>% 
  filter(!taxon %in% c("Homo_sapiens", "Canis_lupus")) %>% # remove these because
  group_by(location, label, type) %>%
  mutate(grp = group_indices()) %>% # this creates a group index for each set of field/pcr replicates
  ungroup() %>%
  select(grp, sample.PCR, taxon, count)

# separate out the grp info
grps <- pool.vegan.df %>%
  select(grp, sample.PCR)

# separate out the df to convert
wide_df <- pool.vegan.df %>%
  select(-grp) %>%
  pivot_wider(id_cols=sample.PCR, names_from=taxon, values_from=count) %>%
  replace(is.na(.), 0) %>% # replace NAs with zeros
  select(-sample.PCR)
  
matrix <- data.matrix(wide_df, rownames.force = T)

convert absolute abundance to relative abundance

# Running NMDS in vegan (metaMDS)
NMS <-
  metaMDS(tax.distmat,
          distance = "bray",
          k = 3,
          maxit = 999, 
          trymax = 500,
          wascores = TRUE)

Run 0 stress 0.1047505 
Run 1 stress 0.1091038 
Run 2 stress 0.1033705 
... New best solution
... Procrustes: rmse 0.01517027  max resid 0.1431347 
Run 3 stress 0.1031345 
... New best solution
... Procrustes: rmse 0.01476544  max resid 0.169073 
Run 4 stress 0.1099409 
Run 5 stress 0.1075893 
Run 6 stress 0.1034186 
... Procrustes: rmse 0.0145414  max resid 0.1562133 
Run 7 stress 0.1030296 
... New best solution
... Procrustes: rmse 0.01489652  max resid 0.1528805 
Run 8 stress 0.1035558 
Run 9 stress 0.1039788 
Run 10 stress 0.1028844 
... New best solution
... Procrustes: rmse 0.01462787  max resid 0.1633978 
Run 11 stress 0.105363 
Run 12 stress 0.1054943 
Run 13 stress 0.1111475 
Run 14 stress 0.1053857 
Run 15 stress 0.1058309 
Run 16 stress 0.1029229 
... Procrustes: rmse 0.01399598  max resid 0.1420114 
Run 17 stress 0.110239 
Run 18 stress 0.1036202 
Run 19 stress 0.1050009 
Run 20 stress 0.103275 
... Procrustes: rmse 0.01339265  max resid 0.164441 
Run 21 stress 0.103434 
Run 22 stress 0.1040243 
Run 23 stress 0.1044319 
Run 24 stress 0.1122533 
Run 25 stress 0.109398 
Run 26 stress 0.1073117 
Run 27 stress 0.1053562 
Run 28 stress 0.1107769 
Run 29 stress 0.1043899 
Run 30 stress 0.1112774 
Run 31 stress 0.1095436 
Run 32 stress 0.1038537 
Run 33 stress 0.1115517 
Run 34 stress 0.1054331 
Run 35 stress 0.1136744 
Run 36 stress 0.1071017 
Run 37 stress 0.106483 
Run 38 stress 0.1033172 
... Procrustes: rmse 0.01208219  max resid 0.1029581 
Run 39 stress 0.1040138 
Run 40 stress 0.1031196 
... Procrustes: rmse 0.01269745  max resid 0.1236347 
Run 41 stress 0.1037689 
Run 42 stress 0.1040558 
Run 43 stress 0.1023241 
... New best solution
... Procrustes: rmse 0.0108321  max resid 0.1406588 
Run 44 stress 0.1075474 
Run 45 stress 0.1042029 
Run 46 stress 0.1075343 
Run 47 stress 0.1036712 
Run 48 stress 0.1044605 
Run 49 stress 0.1093957 
Run 50 stress 0.1047364 
Run 51 stress 0.1071635 
Run 52 stress 0.1040558 
Run 53 stress 0.1069962 
Run 54 stress 0.1078106 
Run 55 stress 0.1035598 
Run 56 stress 0.1032626 
Run 57 stress 0.1042423 
Run 58 stress 0.1104173 
Run 59 stress 0.1068436 
Run 60 stress 0.1037441 
Run 61 stress 0.1066773 
Run 62 stress 0.1064878 
Run 63 stress 0.1139548 
Run 64 stress 0.108836 
Run 65 stress 0.1031774 
Run 66 stress 0.1035821 
Run 67 stress 0.1041209 
Run 68 stress 0.1110943 
Run 69 stress 0.1066009 
Run 70 stress 0.1110221 
Run 71 stress 0.1089743 
Run 72 stress 0.1136639 
Run 73 stress 0.107876 
Run 74 stress 0.1046272 
Run 75 stress 0.1141784 
Run 76 stress 0.1105466 
Run 77 stress 0.1040599 
Run 78 stress 0.1131514 
Run 79 stress 0.1081336 
Run 80 stress 0.1058468 
Run 81 stress 0.1072477 
Run 82 stress 0.1032487 
Run 83 stress 0.1124266 
Run 84 stress 0.1053103 
Run 85 stress 0.1128828 
Run 86 stress 0.1033269 
Run 87 stress 0.1073788 
Run 88 stress 0.1075338 
Run 89 stress 0.1067032 
Run 90 stress 0.1033845 
Run 91 stress 0.10491 
Run 92 stress 0.1053343 
Run 93 stress 0.1038026 
Run 94 stress 0.1039996 
Run 95 stress 0.1057072 
Run 96 stress 0.1043389 
Run 97 stress 0.1031114 
Run 98 stress 0.1085936 
Run 99 stress 0.1035071 
Run 100 stress 0.1085043 
Run 101 stress 0.1043153 
Run 102 stress 0.104713 
Run 103 stress 0.1077827 
Run 104 stress 0.1062981 
Run 105 stress 0.1031222 
Run 106 stress 0.1084776 
Run 107 stress 0.1057907 
Run 108 stress 0.1039365 
Run 109 stress 0.1091723 
Run 110 stress 0.1080755 
Run 111 stress 0.1034593 
Run 112 stress 0.1061215 
Run 113 stress 0.1108576 
Run 114 stress 0.1059668 
Run 115 stress 0.1055788 
Run 116 stress 0.1115548 
Run 117 stress 0.1043042 
Run 118 stress 0.1042351 
Run 119 stress 0.1074185 
Run 120 stress 0.103832 
Run 121 stress 0.1032032 
Run 122 stress 0.1099107 
Run 123 stress 0.1068055 
Run 124 stress 0.1045153 
Run 125 stress 0.1058854 
Run 126 stress 0.1085682 
Run 127 stress 0.1034648 
Run 128 stress 0.1053297 
Run 129 stress 0.1082581 
Run 130 stress 0.103385 
Run 131 stress 0.1038871 
Run 132 stress 0.1031488 
Run 133 stress 0.1099827 
Run 134 stress 0.1115565 
Run 135 stress 0.10442 
Run 136 stress 0.1040222 
Run 137 stress 0.10662 
Run 138 stress 0.1056669 
Run 139 stress 0.1067992 
Run 140 stress 0.106493 
Run 141 stress 0.1083335 
Run 142 stress 0.1067389 
Run 143 stress 0.1070057 
Run 144 stress 0.1079087 
Run 145 stress 0.1059276 
Run 146 stress 0.103159 
Run 147 stress 0.1092317 
Run 148 stress 0.1073626 
Run 149 stress 0.1095654 
Run 150 stress 0.1068984 
Run 151 stress 0.1079596 
Run 152 stress 0.1110705 
Run 153 stress 0.1041155 
Run 154 stress 0.1042449 
Run 155 stress 0.1137419 
Run 156 stress 0.1120904 
Run 157 stress 0.1068519 
Run 158 stress 0.1026228 
... Procrustes: rmse 0.008829727  max resid 0.1013661 
Run 159 stress 0.1100415 
Run 160 stress 0.1044317 
Run 161 stress 0.10397 
Run 162 stress 0.1079679 
Run 163 stress 0.105392 
Run 164 stress 0.1040166 
Run 165 stress 0.1043113 
Run 166 stress 0.1038135 
Run 167 stress 0.1030774 
Run 168 stress 0.1028823 
Run 169 stress 0.1104579 
Run 170 stress 0.1105455 
Run 171 stress 0.1032114 
Run 172 stress 0.1088956 
Run 173 stress 0.1080987 
Run 174 stress 0.1042083 
Run 175 stress 0.1040545 
Run 176 stress 0.1055706 
Run 177 stress 0.1078103 
Run 178 stress 0.1039281 
Run 179 stress 0.104422 
Run 180 stress 0.1038253 
Run 181 stress 0.1050384 
Run 182 stress 0.1060571 
Run 183 stress 0.1089949 
Run 184 stress 0.1052862 
Run 185 stress 0.1040233 
Run 186 stress 0.1131226 
Run 187 stress 0.1100965 
Run 188 stress 0.1065048 
Run 189 stress 0.1041713 
Run 190 stress 0.1123618 
Run 191 stress 0.1041876 
Run 192 stress 0.1067663 
Run 193 stress 0.1029046 
Run 194 stress 0.1080911 
Run 195 stress 0.1086436 
Run 196 stress 0.1109009 
Run 197 stress 0.1078054 
Run 198 stress 0.1092383 
Run 199 stress 0.113815 
Run 200 stress 0.1102316 
Run 201 stress 0.1053111 
Run 202 stress 0.1077149 
Run 203 stress 0.1078718 
Run 204 stress 0.1041129 
Run 205 stress 0.1097915 
Run 206 stress 0.1066643 
Run 207 stress 0.1035781 
Run 208 stress 0.1045359 
Run 209 stress 0.1087377 
Run 210 stress 0.1071641 
Run 211 stress 0.1132297 
Run 212 stress 0.1130048 
Run 213 stress 0.1056566 
Run 214 stress 0.1122072 
Run 215 stress 0.109456 
Run 216 stress 0.1037874 
Run 217 stress 0.1073794 
Run 218 stress 0.1035854 
Run 219 stress 0.1061877 
Run 220 stress 0.108798 
Run 221 stress 0.1062357 
Run 222 stress 0.1040098 
Run 223 stress 0.1132609 
Run 224 stress 0.114126 
Run 225 stress 0.1065758 
Run 226 stress 0.1059224 
Run 227 stress 0.1041753 
Run 228 stress 0.1033392 
Run 229 stress 0.1034376 
Run 230 stress 0.1112893 
Run 231 stress 0.1034152 
Run 232 stress 0.1147076 
Run 233 stress 0.1071721 
Run 234 stress 0.1042944 
Run 235 stress 0.1050444 
Run 236 stress 0.108261 
Run 237 stress 0.1065422 
Run 238 stress 0.1023964 
... Procrustes: rmse 0.01074999  max resid 0.1310516 
Run 239 stress 0.1043985 
Run 240 stress 0.108586 
Run 241 stress 0.1065348 
Run 242 stress 0.1071631 
Run 243 stress 0.110684 
Run 244 stress 0.1083524 
Run 245 stress 0.1053337 
Run 246 stress 0.1078528 
Run 247 stress 0.104499 
Run 248 stress 0.103487 
Run 249 stress 0.1047615 
Run 250 stress 0.1035198 
Run 251 stress 0.103747 
Run 252 stress 0.1058733 
Run 253 stress 0.1034251 
Run 254 stress 0.1067844 
Run 255 stress 0.1086243 
Run 256 stress 0.1045882 
Run 257 stress 0.1063707 
Run 258 stress 0.1030918 
Run 259 stress 0.1040189 
Run 260 stress 0.1082816 
Run 261 stress 0.1095282 
Run 262 stress 0.1059059 
Run 263 stress 0.1084944 
Run 264 stress 0.1046809 
Run 265 stress 0.1096994 
Run 266 stress 0.1093923 
Run 267 stress 0.1028253 
Run 268 stress 0.1149591 
Run 269 stress 0.1047812 
Run 270 stress 0.1078624 
Run 271 stress 0.105012 
Run 272 stress 0.1125052 
Run 273 stress 0.1057068 
Run 274 stress 0.1092531 
Run 275 stress 0.104797 
Run 276 stress 0.1031935 
Run 277 stress 0.1078094 
Run 278 stress 0.1043631 
Run 279 stress 0.1092091 
Run 280 stress 0.1044983 
Run 281 stress 0.1041113 
Run 282 stress 0.1053773 
Run 283 stress 0.1034588 
Run 284 stress 0.1083526 
Run 285 stress 0.1055744 
Run 286 stress 0.1036517 
Run 287 stress 0.1043398 
Run 288 stress 0.1053034 
Run 289 stress 0.1067655 
Run 290 stress 0.105802 
Run 291 stress 0.107568 
Run 292 stress 0.1103305 
Run 293 stress 0.112835 
Run 294 stress 0.1053336 
Run 295 stress 0.1037726 
Run 296 stress 0.1029617 
Run 297 stress 0.1033539 
Run 298 stress 0.1059967 
Run 299 stress 0.1110852 
Run 300 stress 0.1046125 
Run 301 stress 0.1064393 
Run 302 stress 0.1080683 
Run 303 stress 0.1063069 
Run 304 stress 0.1054926 
Run 305 stress 0.1116901 
Run 306 stress 0.103548 
Run 307 stress 0.1081589 
Run 308 stress 0.1125721 
Run 309 stress 0.1097499 
Run 310 stress 0.1039662 
Run 311 stress 0.1147889 
Run 312 stress 0.1046511 
Run 313 stress 0.1049009 
Run 314 stress 0.1086474 
Run 315 stress 0.1072537 
Run 316 stress 0.1032225 
Run 317 stress 0.1032559 
Run 318 stress 0.1037649 
Run 319 stress 0.1054288 
Run 320 stress 0.1104593 
Run 321 stress 0.107653 
Run 322 stress 0.110044 
Run 323 stress 0.1051683 
Run 324 stress 0.1029237 
Run 325 stress 0.1060015 
Run 326 stress 0.1123066 
Run 327 stress 0.1143683 
Run 328 stress 0.1078046 
Run 329 stress 0.1082902 
Run 330 stress 0.1055701 
Run 331 stress 0.1128816 
Run 332 stress 0.1035831 
Run 333 stress 0.1087047 
Run 334 stress 0.1041095 
Run 335 stress 0.1034921 
Run 336 stress 0.1131037 
Run 337 stress 0.1052898 
Run 338 stress 0.1090951 
Run 339 stress 0.1095752 
Run 340 stress 0.1091098 
Run 341 stress 0.1047604 
Run 342 stress 0.1090653 
Run 343 stress 0.1047757 
Run 344 stress 0.1052119 
Run 345 stress 0.1056775 
Run 346 stress 0.1064509 
Run 347 stress 0.1060602 
Run 348 stress 0.1031868 
Run 349 stress 0.1053675 
Run 350 stress 0.1064393 
Run 351 stress 0.1037076 
Run 352 stress 0.1047263 
Run 353 stress 0.1077758 
Run 354 stress 0.1028465 
Run 355 stress 0.1039849 
Run 356 stress 0.1124709 
Run 357 stress 0.1051077 
Run 358 stress 0.1029264 
Run 359 stress 0.1031474 
Run 360 stress 0.1056758 
Run 361 stress 0.1032205 
Run 362 stress 0.1089219 
Run 363 stress 0.1055382 
Run 364 stress 0.1030962 
Run 365 stress 0.1029831 
Run 366 stress 0.1110793 
Run 367 stress 0.1091858 
Run 368 stress 0.1037874 
Run 369 stress 0.1032391 
Run 370 stress 0.1078163 
Run 371 stress 0.1069993 
Run 372 stress 0.1116935 
Run 373 stress 0.1043484 
Run 374 stress 0.1043161 
Run 375 stress 0.1086366 
Run 376 stress 0.1041711 
Run 377 stress 0.1046042 
Run 378 stress 0.1097481 
Run 379 stress 0.1080457 
Run 380 stress 0.1043261 
Run 381 stress 0.103593 
Run 382 stress 0.1047723 
Run 383 stress 0.1050301 
Run 384 stress 0.1081724 
Run 385 stress 0.1137157 
Run 386 stress 0.1048735 
Run 387 stress 0.1060884 
Run 388 stress 0.1059132 
Run 389 stress 0.103554 
Run 390 stress 0.1102933 
Run 391 stress 0.1060102 
Run 392 stress 0.1048778 
Run 393 stress 0.108088 
Run 394 stress 0.1124055 
Run 395 stress 0.1073617 
Run 396 stress 0.1131979 
Run 397 stress 0.1058047 
Run 398 stress 0.1035497 
Run 399 stress 0.1050533 
Run 400 stress 0.1060043 
Run 401 stress 0.1032813 
Run 402 stress 0.103049 
Run 403 stress 0.1047513 
Run 404 stress 0.1063919 
Run 405 stress 0.1062946 
Run 406 stress 0.1099276 
Run 407 stress 0.1050615 
Run 408 stress 0.1040607 
Run 409 stress 0.1038954 
Run 410 stress 0.1067165 
Run 411 stress 0.1093918 
Run 412 stress 0.1041121 
Run 413 stress 0.1063597 
Run 414 stress 0.1060783 
Run 415 stress 0.1045171 
Run 416 stress 0.1041195 
Run 417 stress 0.1105354 
Run 418 stress 0.1082509 
Run 419 stress 0.1095471 
Run 420 stress 0.1064907 
Run 421 stress 0.106217 
Run 422 stress 0.1044849 
Run 423 stress 0.1044891 
Run 424 stress 0.1039686 
Run 425 stress 0.1034413 
Run 426 stress 0.1065345 
Run 427 stress 0.1051907 
Run 428 stress 0.105302 
Run 429 stress 0.103088 
Run 430 stress 0.1078367 
Run 431 stress 0.1081468 
Run 432 stress 0.103194 
Run 433 stress 0.1085562 
Run 434 stress 0.1101137 
Run 435 stress 0.1076278 
Run 436 stress 0.1041889 
Run 437 stress 0.1036967 
Run 438 stress 0.1035539 
Run 439 stress 0.1097133 
Run 440 stress 0.1060898 
Run 441 stress 0.1042543 
Run 442 stress 0.1085263 
Run 443 stress 0.1035459 
Run 444 stress 0.1078768 
Run 445 stress 0.1079541 
Run 446 stress 0.1031527 
Run 447 stress 0.1126149 
Run 448 stress 0.103595 
Run 449 stress 0.1112164 
Run 450 stress 0.106703 
Run 451 stress 0.1051164 
Run 452 stress 0.1058897 
Run 453 stress 0.1037466 
Run 454 stress 0.1048177 
Run 455 stress 0.1124754 
Run 456 stress 0.1081922 
Run 457 stress 0.1066093 
Run 458 stress 0.1041059 
Run 459 stress 0.1039044 
Run 460 stress 0.1035431 
Run 461 stress 0.1072625 
Run 462 stress 0.1080013 
Run 463 stress 0.1083675 
Run 464 stress 0.1041363 
Run 465 stress 0.1043774 
Run 466 stress 0.1028718 
Run 467 stress 0.1057578 
Run 468 stress 0.1043604 
Run 469 stress 0.1041471 
Run 470 stress 0.1079567 
Run 471 stress 0.1085731 
Run 472 stress 0.1084608 
Run 473 stress 0.1071668 
Run 474 stress 0.1082753 
Run 475 stress 0.1040823 
Run 476 stress 0.1052025 
Run 477 stress 0.1109382 
Run 478 stress 0.1045455 
Run 479 stress 0.1064157 
Run 480 stress 0.1027789 
... Procrustes: rmse 0.01328367  max resid 0.1524937 
Run 481 stress 0.1081466 
Run 482 stress 0.1088787 
Run 483 stress 0.1037462 
Run 484 stress 0.1078385 
Run 485 stress 0.1027742 
... Procrustes: rmse 0.01174669  max resid 0.1269963 
Run 486 stress 0.1043833 
Run 487 stress 0.108456 
Run 488 stress 0.1033895 
Run 489 stress 0.1139604 
Run 490 stress 0.1080004 
Run 491 stress 0.1064829 
Run 492 stress 0.1037161 
Run 493 stress 0.1039652 
Run 494 stress 0.1041698 
Run 495 stress 0.1099811 
Run 496 stress 0.1081474 
Run 497 stress 0.1047244 
Run 498 stress 0.1119344 
Run 499 stress 0.1025052 
... Procrustes: rmse 0.01043141  max resid 0.09655726 
Run 500 stress 0.1038412 
*** No convergence -- monoMDS stopping criteria:
   176: stress ratio > sratmax
   324: scale factor of the gradient < sfgrmin

It is common for NMDS analyses to start by running with 2-dimensions (k), but you want to increase the number of dimensions to ensure a minimized stress value. Keep in mind that anything more than 5-dimensions makes it difficult to interpret a 2-dimensional plot.

As a rule of thumb literature has identified the following cut-off values for stress-level:

Higher than 0.2 is poor (risks for false interpretation). 0.1 - 0.2 is fair (some distances can be misleading for interpretation). 0.05 - 0.1 is good (can be confident in inferences from plot). Less than 0.05 is excellent (this can be rare).

A few options: 1) increase k 2) increase maxit

let’s try increasing k

# Running NMDS in vegan (metaMDS)
NMS.k4 <-
  metaMDS(tax.distmat,
          distance = "bray",
          k = 4,
          maxit = 999, 
          trymax = 500,
          wascores = TRUE)

Run 0 stress 0.07386974 
Run 1 stress 0.07491524 
Run 2 stress 0.07437815 
Run 3 stress 0.07403766 
... Procrustes: rmse 0.01514743  max resid 0.1544086 
Run 4 stress 0.07476686 
Run 5 stress 0.07504459 
Run 6 stress 0.07497444 
Run 7 stress 0.0740845 
... Procrustes: rmse 0.01147333  max resid 0.1218068 
Run 8 stress 0.07525555 
Run 9 stress 0.07520172 
Run 10 stress 0.07445861 
Run 11 stress 0.07519227 
Run 12 stress 0.07405646 
... Procrustes: rmse 0.01127032  max resid 0.128445 
Run 13 stress 0.07502137 
Run 14 stress 0.07589977 
Run 15 stress 0.07495102 
Run 16 stress 0.07395249 
... Procrustes: rmse 0.01506798  max resid 0.1277471 
Run 17 stress 0.07414531 
... Procrustes: rmse 0.01108162  max resid 0.1278024 
Run 18 stress 0.07411322 
... Procrustes: rmse 0.01049023  max resid 0.1345969 
Run 19 stress 0.07408317 
... Procrustes: rmse 0.01416774  max resid 0.1311196 
Run 20 stress 0.07425687 
... Procrustes: rmse 0.008242095  max resid 0.1099863 
Run 21 stress 0.07498896 
Run 22 stress 0.07497087 
Run 23 stress 0.07466173 
Run 24 stress 0.0744282 
Run 25 stress 0.07432908 
... Procrustes: rmse 0.01199421  max resid 0.1274044 
Run 26 stress 0.07414727 
... Procrustes: rmse 0.0121427  max resid 0.1257852 
Run 27 stress 0.0738444 
... New best solution
... Procrustes: rmse 0.01148858  max resid 0.1131865 
Run 28 stress 0.07412599 
... Procrustes: rmse 0.01206593  max resid 0.127347 
Run 29 stress 0.07386381 
... Procrustes: rmse 0.01037911  max resid 0.1278972 
Run 30 stress 0.07464 
Run 31 stress 0.07401148 
... Procrustes: rmse 0.01283616  max resid 0.1251787 
Run 32 stress 0.07493285 
Run 33 stress 0.07453871 
Run 34 stress 0.07408212 
... Procrustes: rmse 0.00942815  max resid 0.1120016 
Run 35 stress 0.07491675 
Run 36 stress 0.07408543 
... Procrustes: rmse 0.01112055  max resid 0.1149321 
Run 37 stress 0.0753959 
Run 38 stress 0.0750333 
Run 39 stress 0.07411234 
... Procrustes: rmse 0.0122264  max resid 0.1155168 
Run 40 stress 0.07430884 
... Procrustes: rmse 0.01533867  max resid 0.1551778 
Run 41 stress 0.07526799 
Run 42 stress 0.07397533 
... Procrustes: rmse 0.01238658  max resid 0.1158422 
Run 43 stress 0.07423644 
... Procrustes: rmse 0.01243427  max resid 0.1256216 
Run 44 stress 0.07462486 
Run 45 stress 0.07535054 
Run 46 stress 0.07441149 
Run 47 stress 0.07447003 
Run 48 stress 0.07415679 
... Procrustes: rmse 0.01289608  max resid 0.1260629 
Run 49 stress 0.07449526 
Run 50 stress 0.07419937 
... Procrustes: rmse 0.01292097  max resid 0.1286732 
Run 51 stress 0.07500998 
Run 52 stress 0.07336348 
... New best solution
... Procrustes: rmse 0.01031954  max resid 0.1285916 
Run 53 stress 0.0742318 
Run 54 stress 0.07603917 
Run 55 stress 0.07387452 
Run 56 stress 0.07489372 
Run 57 stress 0.07403113 
Run 58 stress 0.07486089 
Run 59 stress 0.07417056 
Run 60 stress 0.07422359 
Run 61 stress 0.07418035 
Run 62 stress 0.07410374 
Run 63 stress 0.07430432 
Run 64 stress 0.07485444 
Run 65 stress 0.07444278 
Run 66 stress 0.07511817 
Run 67 stress 0.07379722 
... Procrustes: rmse 0.01075407  max resid 0.1286048 
Run 68 stress 0.07528454 
Run 69 stress 0.07419061 
Run 70 stress 0.07373564 
... Procrustes: rmse 0.01101958  max resid 0.1269649 
Run 71 stress 0.07502309 
Run 72 stress 0.07490364 
Run 73 stress 0.07423336 
Run 74 stress 0.07413461 
Run 75 stress 0.07469098 
Run 76 stress 0.0751561 
Run 77 stress 0.07394234 
Run 78 stress 0.0751409 
Run 79 stress 0.07397217 
Run 80 stress 0.07350773 
... Procrustes: rmse 0.008049176  max resid 0.1321206 
Run 81 stress 0.07534628 
Run 82 stress 0.07539287 
Run 83 stress 0.07513191 
Run 84 stress 0.07465574 
Run 85 stress 0.07536908 
Run 86 stress 0.07549339 
Run 87 stress 0.07519677 
Run 88 stress 0.07419013 
Run 89 stress 0.07423067 
Run 90 stress 0.07440885 
Run 91 stress 0.07479592 
Run 92 stress 0.07367861 
... Procrustes: rmse 0.008423985  max resid 0.1316586 
Run 93 stress 0.07577634 
Run 94 stress 0.073537 
... Procrustes: rmse 0.006278642  max resid 0.09073157 
Run 95 stress 0.07456939 
Run 96 stress 0.0742466 
Run 97 stress 0.07426099 
Run 98 stress 0.07533133 
Run 99 stress 0.07511749 
Run 100 stress 0.07384193 
... Procrustes: rmse 0.01201916  max resid 0.1300418 
Run 101 stress 0.0750367 
Run 102 stress 0.0746003 
Run 103 stress 0.07443359 
Run 104 stress 0.07506973 
Run 105 stress 0.07496083 
Run 106 stress 0.0744362 
Run 107 stress 0.07539374 
Run 108 stress 0.07425111 
Run 109 stress 0.07394821 
Run 110 stress 0.07538046 
Run 111 stress 0.07498602 
Run 112 stress 0.07371408 
... Procrustes: rmse 0.009580648  max resid 0.1308994 
Run 113 stress 0.07369853 
... Procrustes: rmse 0.0118808  max resid 0.1286611 
Run 114 stress 0.07520656 
Run 115 stress 0.07487575 
Run 116 stress 0.0742066 
Run 117 stress 0.0741746 
Run 118 stress 0.07447713 
Run 119 stress 0.07550114 
Run 120 stress 0.07531546 
Run 121 stress 0.07412308 
Run 122 stress 0.07405995 
Run 123 stress 0.07536661 
Run 124 stress 0.07518678 
Run 125 stress 0.0750539 
Run 126 stress 0.0742406 
Run 127 stress 0.07527225 
Run 128 stress 0.07410186 
Run 129 stress 0.07486405 
Run 130 stress 0.07510652 
Run 131 stress 0.07417878 
Run 132 stress 0.0743967 
Run 133 stress 0.07428291 
Run 134 stress 0.0734201 
... Procrustes: rmse 0.008337913  max resid 0.09311231 
Run 135 stress 0.07494668 
Run 136 stress 0.07433737 
Run 137 stress 0.07491173 
Run 138 stress 0.07506987 
Run 139 stress 0.07433888 
Run 140 stress 0.07376532 
... Procrustes: rmse 0.01132078  max resid 0.1298581 
Run 141 stress 0.07426083 
Run 142 stress 0.0759338 
Run 143 stress 0.07402262 
Run 144 stress 0.07419204 
Run 145 stress 0.07500366 
Run 146 stress 0.07470428 
Run 147 stress 0.07505205 
Run 148 stress 0.07435951 
Run 149 stress 0.07390386 
Run 150 stress 0.07399403 
Run 151 stress 0.07519878 
Run 152 stress 0.07392592 
Run 153 stress 0.07406892 
Run 154 stress 0.07364193 
... Procrustes: rmse 0.008869131  max resid 0.1300762 
Run 155 stress 0.07389315 
Run 156 stress 0.07494633 
Run 157 stress 0.0741264 
Run 158 stress 0.07520082 
Run 159 stress 0.07526851 
Run 160 stress 0.07400545 
Run 161 stress 0.07511563 
Run 162 stress 0.07518135 
Run 163 stress 0.07443978 
Run 164 stress 0.07461335 
Run 165 stress 0.07505927 
Run 166 stress 0.0746276 
Run 167 stress 0.07392144 
Run 168 stress 0.07497874 
Run 169 stress 0.07538232 
Run 170 stress 0.07505756 
Run 171 stress 0.07522208 
Run 172 stress 0.07469917 
Run 173 stress 0.07500151 
Run 174 stress 0.07389441 
Run 175 stress 0.07385946 
... Procrustes: rmse 0.01400884  max resid 0.1278604 
Run 176 stress 0.07473507 
Run 177 stress 0.07414713 
Run 178 stress 0.07477712 
Run 179 stress 0.07541284 
Run 180 stress 0.07393383 
Run 181 stress 0.07536193 
Run 182 stress 0.07396821 
Run 183 stress 0.0744102 
Run 184 stress 0.07391726 
Run 185 stress 0.0751422 
Run 186 stress 0.07405345 
Run 187 stress 0.07379454 
... Procrustes: rmse 0.01257733  max resid 0.1347662 
Run 188 stress 0.07419108 
Run 189 stress 0.07504814 
Run 190 stress 0.07537465 
Run 191 stress 0.07430705 
Run 192 stress 0.07392313 
Run 193 stress 0.07474243 
Run 194 stress 0.0737347 
... Procrustes: rmse 0.01018345  max resid 0.1163415 
Run 195 stress 0.07393165 
Run 196 stress 0.07441157 
Run 197 stress 0.07399032 
Run 198 stress 0.07382402 
... Procrustes: rmse 0.01144517  max resid 0.1288228 
Run 199 stress 0.0739376 
Run 200 stress 0.07384819 
... Procrustes: rmse 0.01054107  max resid 0.1294516 
Run 201 stress 0.07496302 
Run 202 stress 0.07565388 
Run 203 stress 0.07514652 
Run 204 stress 0.07498086 
Run 205 stress 0.07418292 
Run 206 stress 0.07431398 
Run 207 stress 0.07387517 
Run 208 stress 0.07417817 
Run 209 stress 0.07444683 
Run 210 stress 0.0745562 
Run 211 stress 0.07511673 
Run 212 stress 0.07515831 
Run 213 stress 0.073975 
Run 214 stress 0.07517462 
Run 215 stress 0.07571993 
Run 216 stress 0.07502441 
Run 217 stress 0.07437404 
Run 218 stress 0.07460974 
Run 219 stress 0.07503461 
Run 220 stress 0.07397268 
Run 221 stress 0.07411017 
Run 222 stress 0.07454477 
Run 223 stress 0.07407223 
Run 224 stress 0.0739021 
Run 225 stress 0.07522509 
Run 226 stress 0.0753173 
Run 227 stress 0.07384887 
... Procrustes: rmse 0.0103963  max resid 0.1140664 
Run 228 stress 0.07530704 
Run 229 stress 0.07509589 
Run 230 stress 0.07522791 
Run 231 stress 0.0754115 
Run 232 stress 0.07349118 
... Procrustes: rmse 0.006143577  max resid 0.09491365 
Run 233 stress 0.07404312 
Run 234 stress 0.07492119 
Run 235 stress 0.07512758 
Run 236 stress 0.07402668 
Run 237 stress 0.07520954 
Run 238 stress 0.07491775 
Run 239 stress 0.0755523 
Run 240 stress 0.07419907 
Run 241 stress 0.07542771 
Run 242 stress 0.07423151 
Run 243 stress 0.074433 
Run 244 stress 0.07507299 
Run 245 stress 0.0749457 
Run 246 stress 0.07425704 
Run 247 stress 0.07542338 
Run 248 stress 0.07524835 
Run 249 stress 0.07480125 
Run 250 stress 0.07424095 
Run 251 stress 0.07376875 
... Procrustes: rmse 0.009751902  max resid 0.1187984 
Run 252 stress 0.07505915 
Run 253 stress 0.07394313 
Run 254 stress 0.07429427 
Run 255 stress 0.07536892 
Run 256 stress 0.07511498 
Run 257 stress 0.07420732 
Run 258 stress 0.07398179 
Run 259 stress 0.07429833 
Run 260 stress 0.07547262 
Run 261 stress 0.07428449 
Run 262 stress 0.07429234 
Run 263 stress 0.07444704 
Run 264 stress 0.07448891 
Run 265 stress 0.07415608 
Run 266 stress 0.07527221 
Run 267 stress 0.07408175 
Run 268 stress 0.07509631 
Run 269 stress 0.07516579 
Run 270 stress 0.07454137 
Run 271 stress 0.07541818 
Run 272 stress 0.07514797 
Run 273 stress 0.07396484 
Run 274 stress 0.07521295 
Run 275 stress 0.07497417 
Run 276 stress 0.07469073 
Run 277 stress 0.07373416 
... Procrustes: rmse 0.009084333  max resid 0.1202798 
Run 278 stress 0.07413188 
Run 279 stress 0.07435698 
Run 280 stress 0.07533741 
Run 281 stress 0.07364552 
... Procrustes: rmse 0.01196946  max resid 0.1404214 
Run 282 stress 0.07493169 
Run 283 stress 0.07467212 
Run 284 stress 0.07492725 
Run 285 stress 0.07405564 
Run 286 stress 0.07433395 
Run 287 stress 0.07385769 
... Procrustes: rmse 0.01345113  max resid 0.1259407 
Run 288 stress 0.07371607 
... Procrustes: rmse 0.01199199  max resid 0.1301328 
Run 289 stress 0.07391625 
Run 290 stress 0.07411003 
Run 291 stress 0.07520418 
Run 292 stress 0.07394403 
Run 293 stress 0.07509057 
Run 294 stress 0.0752625 
Run 295 stress 0.07539846 
Run 296 stress 0.07518962 
Run 297 stress 0.07423096 
Run 298 stress 0.07393472 
Run 299 stress 0.07451568 
Run 300 stress 0.07514376 
Run 301 stress 0.07434216 
Run 302 stress 0.07461602 
Run 303 stress 0.07534491 
Run 304 stress 0.07519696 
Run 305 stress 0.07392339 
Run 306 stress 0.07408336 
Run 307 stress 0.07406424 
Run 308 stress 0.07478821 
Run 309 stress 0.07542339 
Run 310 stress 0.07371359 
... Procrustes: rmse 0.01018135  max resid 0.1307369 
Run 311 stress 0.07449296 
Run 312 stress 0.07371978 
... Procrustes: rmse 0.01009987  max resid 0.1289838 
Run 313 stress 0.07386712 
Run 314 stress 0.07395433 
Run 315 stress 0.07472305 
Run 316 stress 0.07481748 
Run 317 stress 0.07489642 
Run 318 stress 0.07549717 
Run 319 stress 0.0750969 
Run 320 stress 0.07406464 
Run 321 stress 0.07406933 
Run 322 stress 0.0744039 
Run 323 stress 0.07542437 
Run 324 stress 0.07416126 
Run 325 stress 0.07527264 
Run 326 stress 0.07370265 
... Procrustes: rmse 0.01107051  max resid 0.1306372 
Run 327 stress 0.0754249 
Run 328 stress 0.07431442 
Run 329 stress 0.07461884 
Run 330 stress 0.07390659 
Run 331 stress 0.07462883 
Run 332 stress 0.07525901 
Run 333 stress 0.07526531 
Run 334 stress 0.07424924 
Run 335 stress 0.0750202 
Run 336 stress 0.07371163 
... Procrustes: rmse 0.01242183  max resid 0.1334272 
Run 337 stress 0.07387081 
Run 338 stress 0.07503239 
Run 339 stress 0.07444035 
Run 340 stress 0.07369797 
... Procrustes: rmse 0.009715762  max resid 0.1364495 
Run 341 stress 0.07490495 
Run 342 stress 0.07452636 
Run 343 stress 0.07514138 
Run 344 stress 0.07466526 
Run 345 stress 0.07460514 
Run 346 stress 0.0751208 
Run 347 stress 0.07405165 
Run 348 stress 0.07468447 
Run 349 stress 0.07545285 
Run 350 stress 0.07399098 
Run 351 stress 0.07506165 
Run 352 stress 0.07545248 
Run 353 stress 0.07495447 
Run 354 stress 0.07387011 
Run 355 stress 0.07458945 
Run 356 stress 0.07422489 
Run 357 stress 0.07517519 
Run 358 stress 0.07422387 
Run 359 stress 0.07493579 
Run 360 stress 0.07424207 
Run 361 stress 0.07466894 
Run 362 stress 0.07519356 
Run 363 stress 0.07375698 
... Procrustes: rmse 0.01079571  max resid 0.1267765 
Run 364 stress 0.07502626 
Run 365 stress 0.07526004 
Run 366 stress 0.07392751 
Run 367 stress 0.07507572 
Run 368 stress 0.07391212 
Run 369 stress 0.07485798 
Run 370 stress 0.07406031 
Run 371 stress 0.07375768 
... Procrustes: rmse 0.01107682  max resid 0.1323648 
Run 372 stress 0.0749523 
Run 373 stress 0.0739031 
Run 374 stress 0.07545146 
Run 375 stress 0.07462628 
Run 376 stress 0.0750187 
Run 377 stress 0.07441731 
Run 378 stress 0.0742416 
Run 379 stress 0.07516328 
Run 380 stress 0.07394378 
Run 381 stress 0.07422517 
Run 382 stress 0.07438329 
Run 383 stress 0.07411848 
Run 384 stress 0.07420582 
Run 385 stress 0.07501456 
Run 386 stress 0.07377241 
... Procrustes: rmse 0.01228797  max resid 0.1278678 
Run 387 stress 0.07410455 
Run 388 stress 0.07519924 
Run 389 stress 0.07393327 
Run 390 stress 0.07504377 
Run 391 stress 0.07561395 
Run 392 stress 0.07353481 
... Procrustes: rmse 0.009013403  max resid 0.1299546 
Run 393 stress 0.07373103 
... Procrustes: rmse 0.01231259  max resid 0.1340134 
Run 394 stress 0.07513421 
Run 395 stress 0.07485496 
Run 396 stress 0.07395286 
Run 397 stress 0.07398958 
Run 398 stress 0.07438072 
Run 399 stress 0.07559357 
Run 400 stress 0.07475406 
Run 401 stress 0.07528713 
Run 402 stress 0.07561433 
Run 403 stress 0.07474959 
Run 404 stress 0.07460552 
Run 405 stress 0.07473478 
Run 406 stress 0.07394682 
Run 407 stress 0.07478529 
Run 408 stress 0.07514026 
Run 409 stress 0.07419717 
Run 410 stress 0.07393309 
Run 411 stress 0.07430879 
Run 412 stress 0.07397659 
Run 413 stress 0.07358595 
... Procrustes: rmse 0.01306045  max resid 0.1258173 
Run 414 stress 0.074179 
Run 415 stress 0.07480772 
Run 416 stress 0.07440817 
Run 417 stress 0.07555891 
Run 418 stress 0.07403137 
Run 419 stress 0.07467252 
Run 420 stress 0.07481558 
Run 421 stress 0.0744667 
Run 422 stress 0.07475577 
Run 423 stress 0.07537888 
Run 424 stress 0.07512247 
Run 425 stress 0.07419997 
Run 426 stress 0.07348721 
... Procrustes: rmse 0.01223881  max resid 0.1276855 
Run 427 stress 0.07408208 
Run 428 stress 0.07399277 
Run 429 stress 0.07443146 
Run 430 stress 0.07525452 
Run 431 stress 0.07369938 
... Procrustes: rmse 0.01272754  max resid 0.1268051 
Run 432 stress 0.07598447 
Run 433 stress 0.07557919 
Run 434 stress 0.07390075 
Run 435 stress 0.07530529 
Run 436 stress 0.07499659 
Run 437 stress 0.07528059 
Run 438 stress 0.07490328 
Run 439 stress 0.07362399 
... Procrustes: rmse 0.01292832  max resid 0.1259931 
Run 440 stress 0.07394942 
Run 441 stress 0.07407381 
Run 442 stress 0.07420546 
Run 443 stress 0.07575533 
Run 444 stress 0.07468921 
Run 445 stress 0.07482684 
Run 446 stress 0.07485404 
Run 447 stress 0.07434552 
Run 448 stress 0.07402306 
Run 449 stress 0.07393731 
Run 450 stress 0.07518266 
Run 451 stress 0.07397705 
Run 452 stress 0.07401406 
Run 453 stress 0.07504258 
Run 454 stress 0.0749351 
Run 455 stress 0.07414586 
Run 456 stress 0.07565358 
Run 457 stress 0.07438571 
Run 458 stress 0.07412992 
Run 459 stress 0.07519675 
Run 460 stress 0.0740214 
Run 461 stress 0.07438382 
Run 462 stress 0.0755578 
Run 463 stress 0.07441952 
Run 464 stress 0.07446079 
Run 465 stress 0.07414533 
Run 466 stress 0.07516437 
Run 467 stress 0.07409466 
Run 468 stress 0.07418167 
Run 469 stress 0.0740646 
Run 470 stress 0.07522238 
Run 471 stress 0.07455339 
Run 472 stress 0.07601397 
Run 473 stress 0.07402227 
Run 474 stress 0.07531606 
Run 475 stress 0.0744596 
Run 476 stress 0.07509474 
Run 477 stress 0.0748124 
Run 478 stress 0.07530191 
Run 479 stress 0.07522585 
Run 480 stress 0.07493452 
Run 481 stress 0.07450819 
Run 482 stress 0.07453053 
Run 483 stress 0.07528691 
Run 484 stress 0.07510803 
Run 485 stress 0.07379753 
... Procrustes: rmse 0.009806388  max resid 0.1298822 
Run 486 stress 0.07544587 
Run 487 stress 0.07397591 
Run 488 stress 0.0748057 
Run 489 stress 0.07371733 
... Procrustes: rmse 0.01116551  max resid 0.1282052 
Run 490 stress 0.07435777 
Run 491 stress 0.07440119 
Run 492 stress 0.07528999 
Run 493 stress 0.07468488 
Run 494 stress 0.07421266 
Run 495 stress 0.07458019 
Run 496 stress 0.07352401 
... Procrustes: rmse 0.006242808  max resid 0.08835656 
Run 497 stress 0.07385832 
... Procrustes: rmse 0.01014896  max resid 0.1302035 
Run 498 stress 0.07417513 
Run 499 stress 0.07425166 
Run 500 stress 0.075162 
*** No convergence -- monoMDS stopping criteria:
   500: scale factor of the gradient < sfgrmin

Still no convergence.

That looks pretty damn good though.

Also not bad.

Let’s go forward with that for the time being.

Plotting the NMDS

# Plotting points in ordination space
plot(NMS, "sites")   # Produces distance 
orditorp(NMS, "sites")   # Gives points labels

BUT, I’m actually less interested in this and more interested in using it with the sample replicates.

Pivoting back to the dataframe

Here is where the uncollapsed df might be more helpful… let me take a quick look to see what Pat did here (and if what I did for Cornell might be useful)